School of Physics and Optoelectronic Engineering, Xidian University, Xi’an 710071, China
† Corresponding author. E-mail:
chengqiouhu@126.com
Project supported by Chinese National Research Fund (Grant No. 9140A02010514DZ01019).
1. IntroductionThe research of numerical methods for optical wave propagation in homogeneous and isotropic media has been an active area of research since the advent of the personal computer.[1] Developments in computational holography have also driven this investigation.[2–6] Although there are several existing numerical methods for different specific situations respectively,[7] researchers are urged to develop faster, more efficient, and more accurate algorithms, in order to satisfy new requirements.[8–11] The angular spectrum method (ASM) is one of the most popular approaches for achieving these objectives. In order to improve the calculation speed, the ASM requires the application of the fast Fourier transform (FFT) and inverse fast Fourier transform (IFFT) instead of traditional numerical integration. However, FFT requires sampling points in an equidistant grid, which may not be satisfied in the case of nonuniform sampling.[12,13]
In this study, an improved ASM algorithm is proposed. The core innovation of this method is that FFT is substituted by trigonometric interpolation (TI). The TI-based ASM (TIASM) is found to be better suited to nonuniform sampling than traditional FFT based ASM (FFTASM). Calculation results also show that TIASM is more amenable to situations where the source field has a simple and strong frequency contrast. For this reason, a billiard + tablecloth is selected as the source field image for the simulation. Irrespective of whether the original phase distribution is a plane wave or a Fresnel zone plate phase distribution, TIASM always produces more accurate results than FFTASM.
2. Traditional angular spectrum methodThe Rayleigh–Sommerfeld integral is a rigorous formula for solving a scalar diffraction problem.[14] When the source plane is parallel to the destination plane as shown in Fig. 1, if the input field is monochromatic, the Rayleigh–Sommerfeld integral is given explicitly by
where
g(
x′,
y′,0) is the distribution of the input field on the source plane,
g(
x,
y,
z) is the distribution of the output field on the destination plane,
z is the distance between the two planes,
r′ = [(
x−
x′)
2 + (
y−
y′)
2 +
z2]
1/2, and
λ is the wavelength.
Instead of integration, equation (1) can be expressed in an alternative concise way based on convolution
where * represents a convolution, and
h(
x,
y,
z) is the propagation kernel expressed as
where
r = [
x2 +
y2 +
z2]
1/2. The equation (
2) can be rewritten using the convolution theorem
where
G and
H represent the Fourier spectrum of functions
g and
h respectively, and
u and
v are the Fourier frequencies along the
x-axis and
y-axis, respectively. Since
u and
v are not the wave number (2
π/
λ) of the optical field, they can be interpreted as the two-dimensional (2-D) signal frequency of the optical field. With the aid of Eq. (
4),
g(
x,
y,
z) can be calculated by
where
f[·] represents a Fourier transform and
f−1[·] represents an inverse Fourier transform. Equation (
5) is the basic equation of the ASM, which is directly derived from the Rayleigh–Sommerfeld integral without any approximation. In numerical simulations,
f[·] and
f−1[·] are replaced by FFT and IFFT respectively, because it is impossible to calculate infinite sampling points from the input field. Owing to this substitution, ASM is often preferred to the Rayleigh–Sommerfeld integral, because computation can be accelerated by using an FFT based algorithm.
3. Trigonometric-interpolation-based angular spectrum methodAlthough the ASM is a popular approach, there is an inherent limitation in this method. For standard ASM, the sampling points are commonly in an equidistant grid in order to facilitate FFT. However, in more general situations, the sampling points may be distributed nonuniformly. For example, a region of the input field may contain local detail with high frequency, while the rest of the field may contain regions of low frequency. It makes FFT not efficient enough due to FFT’s uniform sampling rate through the whole input field.
In this study, a generalized ASM that can deal with nonuniform sampling is proposed. The core idea is that the original FFT process is substituted by a TI. Therefore, this new algorithm is robust to the situation where an input field has variable local frequencies.
In our method, it is assumed that there are N sampling points (xn,yn;gn) in the source field plane, which are nonuniformly distributed. Next, f[g(x,y,0)] is determined via all these nonuniform sampling points. At this point, trigonometric interpolation is implemented. A 2-D function is assumed as p(x,y), which is the sum of trigonometric polynomials
where ±
U and ±
V are the upper and lower bounds of frequencies
u and
v, respectively. The frequency intervals Δ
u and Δ
v are determined by the size of the sampling window as
For the common ASM, the sampling window is actually the source field, which means that both of them have the same size. Both FFT and trigonometric interpolation involve the inherent periodicity of functions in both the Fourier and real space, so there must be errors around the edges of the destination field. In order to eliminate these errors, the area of the sampling window has to be expanded along both the x-and y-axes. For additional areas, extra sampling points must be padded with zeros. In this study, the source field and the sampling window occupy the same center, while the size of the sampling window is twice as long both along the x-and y-axes of the source field. As shown in Fig. 2, we call the additional area around the source field the margin area. After propagation calculations, the output sampling window should be clipped and reduced to the original size again.
Next, we need to calculate coefficients cu,v by using all these nonuniform sampling points. To determine cu,v, every sampling point must satisfy Eq. (6). Therefore, we obtain a set of linear algebraic equations with N individual equations as
Every individual equation above is a trigonometric polynomial and has
M items, where
M satisfy the relation
In Eq. (8), the only remaining unknown is coefficients cu,v. If N > M, there may or may not be a solution of cu,v, depending upon the particular set of sampling points (xn,yn;gn). If N < M, there are infinite sets of solutions of cu,v. If and only if N = M, there is a unique solution of cu,v which satisfies Eq. (8), meaning that the trigonometric polynomials p(x,y) just go through all the sampling points. So we need additional sampling points padded with zeros to fulfill N = M. Furthermore, these additional points must be arranged in the margin area.
At last, a set of solutions of cu,v is obtained, which is equivalent to f[g(x,y,0)]. By Eq. (5), we can calculate g(x,y,z) using the set of cu,v.
4. Calculation results and discussionIn this study, we use an image “billiards” as the amplitude distribution of the source plane as shown in Fig. 3(a). This image with 841 × 841 pixels has a simple background, a tablecloth, so it is suitable for testing our proposed new algorithm. In Fig. 3(b), a Fresnel zone plate is prescribed as the phase distribution, where the pixel value 0 corresponds to a phase value of −π and the pixel value 255 corresponds to a phase value of +π. Its analytic expression is exp [20πi(x2 + y2)/(A/4)2]. The numerical simulation parameters are that the wavelength λ is 1064 nm and the size of the image is 100λ × 100λ. The distance z between the source plane and destination plane is changeable.
In Fig. 4, inside of the two rectangular zones are billiards, which are the main information carriers of the image. So the sampling rate in these two zones is higher than that outside. For the rest of the outside area, the tablecloth, the sampling rate will be reduced by variable scales δ.
When z = 50λ and δ = 5, we calculate the field on the destination plane in Fig. 5. Figures 5(a) and 5(b) are results calculated using numerical integration of Eq. (1) as the correct criteria. To ensure that the errors of the numerical integration are sufficiently small, we make it that one pixel corresponds to one sampling point. Figures 5(c) and 5(d) are results calculated using TIASM. Although TIASM samples fewer points in the tablecloth area than the numerical integration method, both results are very similar, which can be interpreted as that the major information in the image is associated with the billiards and not the tablecloth. Figures 5(e) and 5(f) are results calculated using FFTASM. The sampling rate is five times smaller than Figs. 5(a) and 5(b). From these images, we can find that there is a significant difference between Fig. 5(f) and Fig. 5(b), especially in the center area. We can conclude that the reduction of the sampling rate of the billiards influences the final results. The more points that are sampled in the billiards area, the more accurate the results will be. By adjusting the sampling resolution, we can economize the finite sampling points, and achieve more accurate simulation results.
In Fig. 6, we measure the error of TIASM and FFTASM to Rayleigh–Sommerfeld diffraction which is the accurate criteria, using the signal-to-noise ratio (SNR)[15,16]
where
gRS (
x,
y,
z) is the calculation results by Rayleigh–Sommerfeld integration as accurate criteria.
We can see that TIASM is always better than FFTASM, irrespective of the length of the propagation distance.
In Fig. 7, the original phase distribution of the Fresnel zone plate in Fig. 3(b) is substituted by a plane wave phase distribution, where the pixel value 0 corresponds to a phase value of −π and the pixel value 255 corresponds to the phase value of +π. Its analytic expression is exp(2πix/10λ).
When z = 50λ and δ = 5, we calculate the field on the destination plane in Fig. 8, according to the plane wave phase distribution of Fig. 7. Figures 8(a) and 8(b) are results calculated using a numerical integration of Eq. (1) as the correct criteria. Figures 8(c) and 8(d) are results calculated using TIASM. Figures 8(e) and 8(f) are results calculated using FFTASM. The sampling rate is five times smaller than that of Figs. 8(a) and 8(b). From Figs. 5 and 8, we can conclude that irrespective of whether a plane wave or Fresnel zone plate phase distribution is considered, TIASM always yields more accurate results than FFTASM.
In Fig. 9, we also measure the error of TIASM and FFTASM to Rayleigh–Sommerfeld diffraction when the original phase distribution is a plane wave phase distribution. We can see that TIASM is always better than FFTASM, irrespective of whether a plane wave or Fresnel zone plate phase distribution is considered.
5. ConclusionIn this paper, a new algorithm, TIASM, is proposed and compared with traditional FFTASM. TIASM is especially targeted at nonuniform sampling, which means that the source plane has a variable local sampling resolution. In order to obtain a unique set of solutions for cu,v, we make N = M, so the additional zero sampling points are deployed in the margin area. Therefore, the resulting set, cu,v, is equivalent to a Fourier spectrum for the next calculating process. In the simulation, we use a billiard + tablecloth image. The billiard area represents high local frequency, whereas the tablecloth represents low local frequency information. The results show that TIASM is similar to direct numerical integration of the Rayleigh–Sommerfeld formula. However, traditional FFTASM has significant differences from numerical integration. When the distance is increased, the SNR of TIASM increases at a faster rate in comparison with that of FFTASM. This implies that the advantage of TIASM lies in its suitability for long distance calculation, irrespective of whether the original phase distribution is a plane wave or a Fresnel zone plate distribution. Because TIASM can support variable local resolution, it can capture the local detail of the source field. Therefore, TIASM can increase accuracy in important areas of an image, and reduce the cost of extraneous calculations in trivial areas.
However, there are also two disadvantages for TIASM compared to FFTASM. First, it takes more time to calculate Eq. (8) than traditional FFT. Second, for a different source field, a different artificial area partition is needed, which means that TIASM is not a general approach.
